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Abstract 



Q , The increased sensitivity of future radio telescopes will result in requirements for higher dynamic 

H ' 

range within the image as well as better resolution and immunity to interference. In this paper we propose 



a new matrix formulation of the imaging equation in the cases of non co-planar arrays and polarimetric 
measurements. Then we improve our parametric imaging techniques in terms of resolution and estimation 
accuracy. This is done by enhancing both the MVDR parametric imaging, introducing alternative dirty 



> 

(N 

'^ ■ images and by introducing better power estimates based on least squares, with positive semi-definite 

00 

Tij" ' constraints. We also discuss the use of robust Capon beamforming and semi-definite programming for 

C^ I solving the self-calibration problem. Additionally we provide statistical analysis of the bias of the MVDR 

OO , beamformer for the case of moving array, which serves as a first step in analyzing iterative approaches 



such as CLEAN and the techniques proposed in this paper. Finally we demonstrate a full deconvolution 
process based on the parametric imaging techniques and show its improved resolution and sensitivity 



H ' compared to the CLEAN method. 

Keywords: Radio astronomy, synthesis imaging, parametric imaging, minimum variance, robust beam- 
forming, convex optimization, CLEAN. 

L Introduction 

The future of radio astronomical discoveries depends on achieving better resolution and sensitivity 
while maintaining immunity to terrestrial interference which is rapidly growing. To achieve the improved 
sensitivity and higher resolution new instruments are being designed. The Square Kilometer Array (SKA) 
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II m and the Low Frequency Array (LOFAR) [2] are two of these advanced instruments. Achieving 
higher sensitivity to observe faint objects results in high dynamic range requirements within the image, 
where strong sources can affect the imaging of the very weak sources. On the other hand, Moore's law |[3l 
together with recent advances in optimization theory |4| open the way to the application of more advanced 
computational techniques. In contrast to hardware implementation, these image formation algorithms, that 
are implemented in software can benefit from the continuing increase in computational power, even after 
the antennas and the correlators will be built. In this paper we extend the parametric deconvolution 
approach of ||5l to obtain better power estimation accuracy and higher robustness to interference and 
modeling uncertainty. The algorithms presented here can also be used in conjunction with real-time 
interference mitigation techniques as described in Q and m. 

We briefly describe the current status of radio astronomical imaging techniques. For a more extensive 
overview the reader is referred to f/l, ||8l or ||9l. A good historic perspective can be found in lilOJ . whereas 
lOn provides a very recent perspective. 

The principle of radio interferometry has been used in radio astronomy since 1946 when Ryle and 
Vonberg constructed a radio interferometer using dipole antenna arrays [|12il . During the 1950's several 
radio interferometers which use the synthetic aperture created by movable antennas have been constructed. 
In 1962 the principle of aperture synthesis using earth rotation has been proposed |[T3l . The basic idea 
is to exploit the rotation of the earth to obtain denser coverage of the visibility domain (spatial Fourier 
domain). The first instrument to use this principle was the five kilometer Cambridge radio telescope. 
During the 1970's new instruments with large aperture have been constructed. Among these we find the 
Westerbork Synthesis Radio Telescope (WSRT) in the Netherlands and the Very Large Array (VLA) in 
the USA. Recently, the Giant Microwave Telescopes (GMRT) has been constructed in India and the Allen 
Telescope Array (ATA) in the US. Even these instruments subsample the Fourier domain, so that unique 
reconstruction is not possible without some further processing known as deconvolution. The deconvolution 
process uses some a-priori knowledge about the image to remove the effect of "dirty beam" side-lobes. 

Two principles dominate the astronomical imaging deconvolution. The first method was proposed by 
Hogbom 1141 and is known as CLEAN. The CLEAN method is basically a sequential Least-Squares 
(LS) fitting procedure in which the brightest source location and power are estimated. The response of 
this source is removed from the image and then the process continues to find the next brightest source, 
until the residual image is noise-Uke. During the years it has been partially analyzed ifTSl . |[T6l and ifTTl . 

'For information on the SKA project tlie reader is referred to littp://www.skatelescope.org. 
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However full analysis of the method is still lacking due to its iterative nature. The CLEAN algorithm 
has many recent flavors which are capable of faster performance e.g., the Clark version lITSl and the 
Cotton-Schwab ||T9l. In these versions the dirty image is recomputed only after several point sources 
have been estimated. Furthermore the sources are subtracted from the ungridded visibility. This results in 
better suppression of the sources. Interestingly it is well known that when the noise model is non-white 
(as in radar clutter modeled by ARMA processes in SAR applications) algorithms such as RELAX EOl 
outperform the CLEAN algorithm. 

A second approach proposed by Jaynes f2V\ is maximum entropy deconvolution (MEM). The basic 
idea behind MEM is the following. Not all images which are consistent with the measured data and the 
noise distribution satisfy the positivity demand, i.e., the sky brightness is a positive function. Consider 
only those that satisfy the positivity demand. From these select the one that is most likely to have been 
created randomly. This idea has also been proposed in fT2\ and applied to radio astronomical imaging in 
ll23l . Other approaches based on the differential entropy have also been proposed ll24l . |[25l . An extensive 
collection of papers discussing the various methods and aspects of maximum entropy can be found in 
the various papers in |26|. Briggs [|27i proposed a non-negative least squares approach (NNLS) which 
eliminates the need for iterative processing. However, the computational complexity is very large. 

In this paper we use a reformulation of the image formation problem as a parameter estimation problem 
using a set of covariance matrices, measured at the various observation epochs [5]. This yields a model 
where the array response is time varying. Previous research on time varying arrays and their application 
to direction-of arrival (DOA) estimation includes |[28l . |[29l and ||30l . In iQ we proposed a simplified 
ML estimator. Lanterman |[3T| developed a full EM algorithm for implementing the MLE proposed in 
|[5l . The algorithm performs quite well although it is quite complex compared to the solutions described 
in this paper. 

The above mentioned algorithms assume perfect knowledge of the instrumental response (point spread 
function). Due to various internal and external effects this assumption holds only approximately. One way 
to overcome this problem is the use of calibrating sources. An unresolved source with known parameters 
is measured, and by relating the model errors to the array elements a set of calibration equations is 
solved. A much more appealing solution is to try to improve the fitting between the data and the sky 
model by adjusting the calibration parameters. Another possibility |[32ll is to use the redundant structure 
of the array to solve for the calibration parameters (this is possible only for some arrays which have 
redundant baselines, such as the WSRT). A good overview of the various techniques is given in |[33l . 

In this paper we extend the above methods in several directions. First, we extend the parametric formu- 
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lation of 121 to the non-coplanar array case and to polarimetric imaging. Then we propose relatively "low" 
complexity approaches for the image deconvolution problem, based on Minimum Variance Distortionless 
Response (MVDR) and its robust extensions. We call these extensions Least-Squares Minimum Variance 
Imaging (LS-MVI). We provide a new type of dirty image that has isotropic noise response, something 
desirable in imaging applications. This is done by generalizing the work of Borgiotti and Kaplan 1341 
to the moving array case. We discuss acceleration techniques, related to the Clark [1I81I and the Cotton- 
Schwab II19II approaches to CLEAN. However, in contrast to the classical CLEAN case, the accelerated 
algorithm involves semi-definite programming of low order, ensuring that the covariance matrices remain 
positive semi-definite after the subtraction. We provide analytic expressions for the asymptotic bias of 
the MVDR based imaging. We also relate the classical self-calibration technique to a novel extension 
of the robust Capon beamformer to the moving array case, showing that self-calibration can be cast in 
terms of semi-definite programming. This has advantage over previous approach to self-calibration since 
the covariance structure maintains its positive definite structure. 

We also demonstrate full parametric deconvolution process based on the proposed technique and 
compare CLEAN and LS-MVI on simulated images. We will show in simulations that LS-MVI, indeed, 
significantly outperforms the classical CLEAN algorithm over a wide range of parameters, in terms of 
better resolution, higher sensitivity and dynamic range. 

The structure of the paper is as follows. In section JI] we describe the astronomical measurement 
equation. The measurement equation is subsequently rephrased in a more convenient matrix formulation 
both for non co-planar and co-planar arrays. It is then extended with the effect of noise and unknown 
calibration parameters. We also discuss extension to polarimetric measurements. In section |lll] we discuss 
the new Least Squares Minimum Variance Imaging and extensions to self calibration using robust Capon 
techniques. In section |IV] we present bias analysis of the MVDR imaging technique. 

In |V] we describe several computer simulations demonstrating the gain in parametric deconvolution 
relative to the CLEAN in terms of resolution, sensitivity and capability to model extended structures. We 
also provide example of the tightness of the statistical bias analysis of the MVDR DOA estimation with 
a moving array. We end up with some conclusions. 

II. Astronomical measurement equations 

In this section we describe a simplified mathematical model for the astronomical measurement and 
imaging process. Our discussion follows the introduction in l9l and its matrix extension in [5|. We 
extend the matrix formulation of l5l| to non co-planar arrays and polarimetric measurements. This allows 
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us to obtain a uniform description of various astronomical imaging operations such as deconvolution and 
self-calibration. 

A. The interferometric measurement equation 

The waves received from the celestial sphere may be considered as spatially incoherent wideband 
random noise. They may be polarized and can contain spectral absorption or emission lines. Rather than 
considering the emitted electric field at a location on the celestial sphere, astronomers try to recover the 
intensity (or brightness) //(s) in the direction of unit-length vectors s, where / is a specific frequency. 
Let Ef{r) be the received celestial electric field at a location r on earth and ^(s) be the antenna gain 
towards direction s. The measured correlation of the electric fields between two identical sensors i and 
j with locations rj and r^ is called visibility and is (approximately) given by [9o 

Vf{r,,r,):=E{Ef{r,)Ej(:^} = [ A{s)If{s)e-^-^f '^^""^-'^^'^ d^ . (1) 

J sky 

E{ • } is the mathematical expectation operator, the superscript ^ denotes the transpose of a vector, and 

overbar denotes the complex conjugate. Note that it is only dependent on the oriented distance r^ — r^ 

between the two antennas; this vector is called a baseline. 

For simplification, we may sometimes assume that the astronomical sky is a collection of d discrete 

point sources (maybe unresolved). This gives 

d 

If{s) = Y,If{si)6{s-si), 
1=1 
where s/ is the coordinate of the /'th source, and thus 

d 
Vf{r„vj) = J;/;(s0e-2-•'^^^('■'-'■^■)/^ (2) 

1=1 
For earth rotation synthesis arrays, the following coordinate system is often used. We assume an array 

with antennas that have a small field of view and that track a reference source direction sq in the sky. 

Other locations in the field of view can be written as s = sq + cr , sq -L cr (valid for small a) and a 

natural coordinate system is sq = [0, 0, if'^jCr = [£, m, 0]^. Similarly, the receiver baselines can be 

parameterized as rj — rj = X[u, v, w]'^ , -^ = 7 • The measurement equation in {u, v, w) coordinates 

thus becomes (see Perley 121 chapter 19) 

J — oo J—oo V -L — t — Tn 
^To simplify notation, from this point we do not include the directional response of the elements of the radio telescope. 
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Assuming phase tracking is performed to compensate for the geometric delay we obtain 



OO fOO 



V'{U,V,W)= I I /(/,m) ^-2nj{ul+vm+wn)^j^^ (4) 

J~oo J~oo V 1 — /2 — m^ 
where n = a/I — P — m?. Assuming that the sky is composed of discrete set of sources we obtain 

V'{U, V,W)=Y^ ^(^^'""^) ^~2n j{ul,+vm,+wn,) (5^ 

e=i \ I -I'j -mj 



In certain conditions, such as East-West array or when the field of view is limited by the antenna primary 
beam to sufficiently small angular region, the third term in the exponential can be neglected and the term 



I — Ij — rnj is approximately 1. The measurement equation ^ becomes: 

Vf{u,v) = ff If{l,m)e~^''^'^''^+'"^Uldm. (6) 



which has the form of a Fourier transformation. 

The function Vf{u,v) is sampled at various coordinates {u,v) by first of all taking all possible sensor 
pairs i,j or baselines fj — r^, and second by realizing that the sensor locations r^, tj are actually time- 
varying since the earth rotates. Given a sufficient number of samples in the {u, v) domain, the relation 
can be inverted to obtain an image (the 'map'). Direct Fourier inversion of the visibility data suffers 
from severe aliasing and is termed the dirty image. To overcome the aliasing, deconvolution algorithms 
are required. 

B. Matrix formulation of the measurement equation 

To allow a parametric imaging formulation of the discrete source model, we can now formulate our 
measurement equations in terms of matrices following Q. Our formulation is more general then lU, 
extending it to the non co-planar array case. 

Let ro(tA;) be an arbitrary and time- varying reference point, typically at one of the elements of the 
array, and let us take the {u,v,w) coordinates of the other telescopes with respect to this reference, 
rj(t) — ro(t) = A[ujo(i), Vio{t), mjo(t)] , i = !,••• ,P- Similarly to |5] in the coplanar case, 

equation ^ can then be written in matrix form as 

Rfe = AkBAl, (7) 

where 

Rfc = K{tk), Ak = [afc(/i,mi),... ,a.k{ld,md)] , 
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a.kil,m) 



- 27rj(uio (tfc )i+t)io (tfc )m+tOio {tu )n) 



-2-Kj(upa {tk)l+Vpo (tk)m+Wpa (tk)n) 



(8) 



and 



B = diag 



I{li,mi] 



I{ld,nid) 



yr^ 



mt 



1 



/2 



mi 



1 — ij — mj is close to unity and can be neglected. 



When the field of view is limited the factor 
The vector function afc(/,m) is called the array response vector in array signal processing. It describes 
the response of the telescope array to a source in the direction {l,m). As usual, the array response is 
frequency dependent. The response is also slowly time-varying due to earth rotation. It is assumed that 
the function as shown here is completely known. However, typically the array response is not perfectly 
known. Each antenna may have a different complex receiver gain, 7j(t), dependent on cable losses, 
amplifier gains, and (slowly) time varying. For LOFAR type arrays the calibration is also space varying 
because of atmospheric conditions. To simplify the exposition we only treat calibration parameters that 
are spatially invariant assuming that the field of view is contained in a single isoplanatic patch. The case 
of LOFAR type array with space varying calibration parameters is complicated, and will be treated in 
a subsequent paper. The measurements are also containing additive system noise. Typically this noise is 
zero mean, independent among the antennas (thus spatially white). After noise calibration and scaling of 
the measurements we can also assume that it has a covariance matrix that is a multiple of the identity 
matrix, a^I, where a'^ is the noise power on a single antenna. Usually the receiver noise is assumed to 
be Gaussian. The resulting model of the received covariance matrix then becomes 



^k 



TkAkBAlTl + a^I 



where 



diag[7i,fc, 



I Tp,fcJ 



(9) 



(10) 



and 7j fc = ji{tk) is the calibration parameter for antenna i at time epoch tk- To allow proper calibration it 
is assumed that the time varying calibration parameters vary slowly so calibration equations can combine 
multiple epochs. 



C. The astronomical measurement equation for polarized sources 

Next, we extend our matrix formulation to the case of polarimetric measurement of polarized sources. 
Our notation follows the conventions of Hamaker et. al |[35l - ll36l . For concreteness, we concentrate on 
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circular feeds. We perform the deconvolution of the cross-polarization parameters and then recover the 
Stokes parameters using the MuUer matrices for each source. This enables a straightforward extension of 
the matrix formulation to polarized sources. First we restrict our attention to purely orthogonal feeds, and 
then introduce the Jones matrices of each antenna into the model. To that end assume that each antenna 
has two orthogonal circular feeds. As usual we restrict our discussion to the quasi-monochromatic case. 
The response of the array at epoch k = 1, ...,K towards direction s can now be decomposed as 



[aL,fc(s), aR,fc(s)] = afc(s) h 



(11) 



where I2 is the 2x2 identity matrix. This implies that for any directions si, S2 the vectors a/, ^(si), a^ ^(82) 
are orthogonal. Assume that the sky is composed of d point sources (to replace the integral by a finite 
sum). The extended measurement equation for ideal feeds is now described by 



Rt 



A«B(A«)^ + a2l 



where 



A^ = [aL,fe(si),a^_fc(si),...,ai^fc(sd)a/j_fc(srf)] 



is the array response matrix, 



B = diag{Ei,...,Erf} 



(12) 



(13) 



(14) 



hjf 



is a block diagonal matrix with 2x2 sub-blocks, E^, on the diagonal consisting of the source coherency 
matrices 

Similarly to the calibrated model for scalar imaging, we can introduce feed calibration matrices V^ for 
non-ideal (non-orthogonal) feeds, using the Jones matrices: 



(15) 



r = diag{j\...,JP}, 



(16) 



where the 2x2 sub-blocks are the corresponding Jones matrices representing the polarization leakage and 
the different gains and phases of the feeds: 






•^21 ^^22 

The calibrated measurement equation in matrix form now becomes 

Rfc = r,A®B(A®)^rf + ^2j_ 



(17) 



(18) 
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This model is very similar to the scalar measurement equation ^, except that the source matrix is now 
block diagonal. This will have an effect on the parametric imaging algorithms. Note that it is easy to 
see that the 2x2 sub-blocks of the matrix R^ agree with the formulation of Hamaker et. al ||36l . After 
deconvolution of the model (IT2l ) we can extract the Stokes parameters of each source using the MuUer 
matrices S 

he = Sef , (19) 

where b^ = vec (E^) and 

e^ = [Ie,Qi,Ue,Vef. (20) 

III. The Least-Squares Minimum-Variance Imaging 

The idea of using Direction-of-Arrival (DOA) estimation techniques for imaging was first introduced 
in Q. In that paper it was suggested that the imaging process will be based on MVDR (Minimum 
Variance Distortionless Response) DOA estimates |[37l . using classical dirty image intensities, similar to 
the CLEAN method |[T4l . Therefore, the algorithm improves the source location estimates, but the power 
estimates are still inaccurate similarly to the CLEAN algorithm. In this paper we improve the algorithm 
of in in several directions: First, we introduce a new type of dirty image that has better properties than 
the standard MVDR dirty image. We term this new dirty image the Adaptive Angular Response (AAR) 
dirty image, since it extends the technique of |[34ll to the moving array case. The main advantage is the 
isotropic noise response, similar to that of the classical Fourier beamformer. Then, we improve the power 
estimate using non-negative least squares estimator for the power and imposing a positive semi-definite 
constraint on the residual covariance matrices. Note that in contrast to the global NNLS of Briggs ||27l 
we only solve linear NNLS which has a closed form solution. Furthermore, the dimensionality of the 
specific NNLS used in this paper is very low, operating on a single (or few) sources at a time. The two 
previous papers on parametric imaging Q, ll38l have only demonstrated dirty images based on MVDR 
and the Robust MVDR estimation. In this paper we demonstrate the benefits of the parametric approach 
after a full deconvolution. The proposed algorithm results in a robust technique capable of much better 
resolution compared to the CLEAN, since it inherently suppresses interference from other sources within 
the image, by using a data dependent beam in order to estimate the locations. Finally, we will show how 
to speed up the proposed algorithm by estimating multiple sources at each round. 

The MVDR based imaging can be generalized to Robust Capon Beamforming (RCB) 1311, HH, l40l . 
EH, l42l . l43l . which is also expected to provide robustness to array manifold errors. A first stage of 
experimental verification of using the RCB has been carried out in l38l . where dirty images based on the 
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RCB showed superior interference immunity compared to the ordinary dirty images. We also demonstrate 
how to combine the LS-MVI with the self-caUbration. However, this requires a novel extension of the 
robust Capon algorithm to the case of moving array using semi-definite programming. We begin with 
discussion of the MVDR and the robust MVDR algorithms for a moving array and then discuss power 
estimation. Then we provide a complete description of the algorithm and its accelerated versions. We 
end up with extending the robust Capon beamformer to the moving array case and demonstrate how it 
can be used for self calibration. 

A. Two dimensional minimum variance estimation with a moving array 

The Minimum Variance Distortionless Response method ( MVDR ) |[37l was one of the first super- 
resolution techniques for direction-of-arrival estimation. Compared to classical (Fourier) based beam- 
forming it provides better separation of closely spaced sources and robustness to strong interference. The 
MVDR estimator is obtained by solving the following optimization problem p?^: 
Let x(tfc) denote the antennas' output signal at time t^ (we assume that DC has been removed so that 
x(tfc) is a zero mean random signal). We apply a weight vector w^ to x(ifc). Hence, the variance of the 
antennas' output, derived from pointing towards direction s, is given by: 

yar[wf (s)x(tfc)] = wf (s)RfcWfc(s) . 
The total output power is then given by 

I'D{^)-=Y.^k{^)^k^k{^)- (21) 

k 

Note that for the classical dirty image Wfc(s) = afc(s) iQ and the dirty image is given by 

/B(s):=^a^(s)Rfcafc(s). (22) 

k 

This amounts to using fixed Fourier basis vectors for beamforming, independently of the data. As in [5], 
we can replace R by R — cr^I, where cr^ is the noise power. This power downdating can, alternatively, 
be incorporated into the positive semidefinite constraints. A much better approach to the imaging process 
would be to minimize the interference subject to transferring the desired direction unchanged. This is 
equivalent to working with data dependent beamformer. Using the data to form the beam provides much 
better interference suppression, compared to the fixed Fourier basis. The variance of the array output 
consists of the response to sources at many directions. We require that at each time instance k the total 
output power will be minimized, subject to the constraint that the output response towards direction s 
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will be fixed. This is equivalent to requiring that the contributions of the sidelobes of other sources will 
be minimized. Using equation ^ and assuming calibrated array we obtain 

Rfe = AfcBAl^ + a^I . 

In order to obtain an estimate of the power originating from direction s without interference, we require 
that for all k 

wf (s)afc(s) = 1 (23) 

and then minimize the overall output power. This can be reformulated as the following constrained 
problem (For simplicity we denote wa,(s) by Wj.): 

Wfc(s) = argmin w^Rjtviffc subject to w^afc(s) = l. (24) 

Wfc 

The above problem can be solved using Lagrange multipliers. The solution is given by: 

Wfc(s) = ^fc(s)R^^afc(s), where /3fc(s) = . (25) 

afe(s)Rfc afc(s) 

Substituting Wfc(s) in Eq. (|2T| | we obtain 



K 



I'^(s) = y J- . (26) 

Definition 3.1: I'£)(s) formulated in Eq. (l26l) is the MVDR dirty image. 

Finally, we describe a second variant of the MVDR dirty image, which combines the adapted angular 

response (AAR) of 041 and the approach of fl3^ll to moving arrays. Borgiotti and Kaplan 1341 proposed 

to use an MVDR type of estimator w(s) = ;uR~^a(s) but in order to obtain isotropic behavior of the 

noise to add a constraint ||vif(s)|p = 1 for all s. This results in the following spatial spectrum estimator 

a(s)^R-ia(s) 

^ ^ , ^. (27) 

a(s)^R-2a(s) 

Rieken and Fuhrmann ll45l suggested to assume that the measurements at each epoch (or array in their 

formulation) are uncorrelated, resulting in block diagonal covariance matrix 

R = diag{Ri,...,Rfc}. (28) 

This assumption holds for the radio-astronomical case, since the received signal can be assumed inde- 
pendent over different time epochs. Let 

w = ^R^^a, (29) 
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where, 



a 



[ai(s) ,...,aii-(s) 



/;w = f:r " ' "^. : ! ■ oo) 



Substituting (1281 ) and ( |29l ) into (I27l ). we obtain that the dirty image is given by 

Ef=ia,(s)^R^^a,(s) 

This new spatial power spectrum estimator for the moving array case has isotropic white noise response, 
and optimal suppression of interference under the isotropic white noise requirements. Furthermore, the 
averaging over the time epochs k = 1,...,K at both numerator and denominator results in smoother 
behavior of the dirty image. 

Definition 3.2: I'^(s) formulated in Eq. (l30l ) is the AAR dirty image. 
The AAR dirty image (l30l ) yields better power estimates and is more robust to strong noise, but has 
somewhat higher computational complexity. 

B. The LS-MVI Algorithm 

In this section we present a new algorithm - The Least Squares Minimum Variance Imaging (LS- 
MVI). The algorithm is an iterative one, similar to the CLEAN method. We assume that every cosmic 
source in the sky brings the MVDR dirty image ( [261 ) or the AAR dirty image (|30l ) to its maximum at its 
direction. Thus, in each iteration we find the brightest point in the MVDR dirty image. Then we estimate 
its intensity, using Least Squares, as described below. We subtract part of the source's contribution to the 
correlation matrices. A new MVDR dirty image is then calculated, using the new correlation matrices. 
We continue the iterations until a certain stopping rule (typically defined by x^ test for the residual dirty 
image) is met. The final image is composed of the locations and intensities we have found during these 
iterations, convolved with a synthesized beam, usually an ideal Gaussian beam. Note that the location 
estimator is not limited to the grid of the dirty image. Since the model is continuous either interpolation 
or local optimization can be used to find the best location, independent of the grid. Typically, quadratic 
interpolation around the maximum of the grid suffices. This solves the dynamic range problem pointed 
out by B6l . In lH the locations of the sources were estimated by using the MVDR dirty image /^(s). 
However, the intensities were estimated using the conventional dirty image /d(s). As explained before 
the new AAR dirty image improves the location estimation over the MVDR dirty image. Furthermore, we 
go even beyond the location estimation, by improving the estimate of the sources intensity. We suggest 
to estimate the intensities by using Least Squares. Recall from (O that 
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In each iteration we find a location sq that brings the MVDR dirty image to its maximum at its direction. 
Thus, defining a = /(so), we solve the following problem: 

a = ar:gmm^^yJ2k=i \\^k-aa.k{so)EL^{so)\\p, s.t.a>0. (31) 

where || • \\f is the Frobenius norm. Problem (|3TI) can be reformulated as follows: 

a = argmin ||x — ah||2 , s.t. a>0, (32) 

{a} 

where 

x=Kc(R«)^...i;ec(RW)T 

is a vector of size (Kp'^xl), in which K is number of time instances and p is the number of antennas, 
the index ^^"^ indicates the specific time instance, and 

'a(i)(so) a(i)(so))^ , . . . , (a(^)(so) ® a(^)(so)^ " 

The problem can be solved using the Karush-Kuhn-Tucker Conditions for constrained optimization [4]. 

The result is 

fh^x 1 
« = max|^-^,0| . (33) 

When positivity constraint is not applicable, such as in polarimetric imaging, one similarly obtains a = 
jpT^. The MVDR dirty image /^(s) can be replaced with the AAR dirty image (l30l ) yielding better 
performance. The superiority of the LS-MVI over the CLEAN method is shown in simulated examples 
in section IV-AI A possible improvement of the method above is by using semi-definite constraints. We 
can add the following constraints on a in ( [32l ): 

Rfc-cj2l-aafe(s)afc(s)^ ^ k = l,...,K (34) 



where A ^ B means that B — A is positive semi definite. Since the solution for ([32]) provides an upper 
bound on a and is a lower bound, a simple bi-section can provide this optimal value. The LS-MVI 
algorithm is described in Table IH 

C. Accelerating the algorithm 

The proposed algorithm has higher computational complexity compared to the CLEAN since it cannot 
utilize the two-dimensional EFT algorithm to generate the MVDR dirty image. It should be noted that 
naive implementation of the techniques is much more complex than the CLEAN approach. It is still an 
interesting research problem to reduce the complexity of the proposed parametric approach. There are 
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TABLE I 
The LS-MVI algorithm 



n = 0. 








Set Rf^ = Rfc. 








Set the loop gain 7 (typically 7 = 0.1). 








Calculate I^ according to Eq. J30|l. 








last iteration flag=false. 








While not last iteration flag 








s(") = argmax{,j (|I^|). 








Calculate a„ accordine to Eqs. ll33k and ll32k 






Improve the estimate of a„ using J34|(. 








Ri"+^' = Ri") -7ana4s("')afe(s(" 


)", 


k-- 


= 1,...,K. 


Calculate I'^ using Ri"+^'. 








n=n+l. 








Compute last iteration flag (Typically usir 


'gx' 


est 


on residual image. 


or using MDL). 








End. 








I = I^n7 «« Bsynthil -ln,m-m„). 









however several techniques that can accelerate the LS-MVI algorithm. We would like to mention two 
approaches. First a clever update of the MVDR dirty image is possible, with complexity KN'^p^, where 
p is the number of antennas in the array, A^^ is the image size and K is the number of observation 
epochs, by using low rank updates. Moreover, since the MVDR estimator better suppresses the effect of 
other sources, we can estimate simultaneously several point sources with directional vectors si,...,Srf, 
where d is a small number, by choosing the d strongest points in the image. This is similar in spirit to 
the small cycle in lITSi and |[T9l . To that end define 



i«(s,)55a«(s,) 



T 



iW(s.)®aW( 



-\T 



and let H = [hi , . . . , h^] . We can estimate the powers of the sources by solving the problem 



a = argmin ||x — Hq;| 
{a.} 



s.t. Qj > ,« = 1, ..., d 



(35) 



(36) 



where a = [ai^. . . ,ad] are the sources powers. Similarly to the single source case, the solution 
is given by solving the LS problem with positivity constraints. This problem is a special case of 
quadratic programming problems, and therefore can be solved efficiently. In contrast to [27 1, where all the 
measurement constraints are put into a large constrained LS problem, our problem has low dimensional 
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positivity constraints. The solution is given by 

ai = max | ('(H^H) "^ H^x) , o} . 

Therefore, we solve an unconstrained LS solution and set the negative terms to 0. Solving for d sources 
simultaneously reduces the complexity by a factor of d. In this case we can also impose the extra 



conditions 

,(m) _ 

i=l 



^t^ - Yl «^afc(si)afc(s,)^ ^0, k = l, .., K. (37) 



These constraints prevent the matrices R^™ from becoming non positive semi-definite and prevent 
overestimation of the power. This problem is equivalent to a special case of convex programming called 
semi-definite programming in the variables ct = [qi, ..., a^], t 

mill t 
s.t. 

I (x - Ha) 



(x - Uaf t 



>- 



(38) 
^ti"iafc(si)a|^(si) ^ Rfe k = l,...,K 

< t 

After a = [ai, ...,ad] is estimated, the covariance matrices R^ are updated, by subtracting the contri- 
bution of the estimated sources from the covariance matrices (possibly using a predefined loop gain) 

d 

^(n+1) ^ ^(n) _^j2 a,a,(s,)af (s,) (39) 

i=l 

Such problems can be solved very efficiently using convex optimization techniques 101, especially, when 
the number of sources estimated simultaneously is small. Finally, we mention that the MVDR dirty image 
can be interpolated in the intermediate steps, resulting in significant reduction in the complexity of the 
deconvolution. 

D. Self calibration and robust MVDR with a moving array 

We now turn to the case where the array response is not completely known, but we have some statistical 
knowledge of the error, e.g., we know the covariance matrix of the array response error at each epoch. 
Typically this covariance will be time invariant or will have slow temporal variation. In this case we 
extend the robust dirty image as described in [381 . into the moving array case. This generalization is 
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new, and has not been previously dealt with in the signal processing literature. Since the positive definite 
constraint on the residual covariance matrices is important in our application, we decided to extend the 
robust Capon estimator of |39'|. To that end assume that at each epoch we have an uncertainty ellipsoid 
describing the uncertainty of the array response (as well as unknown atmospheric attenuation). This is 
described by 



\H 



(afc(s) - afc(s))^ Cfc (afc(s) - afc(s)) < 1 



(40) 



where afc(s) is the nominal value of the array response towards the point s. Generalizing our MVDR 
with moving array we would like to solve the following problem: 

[p,ai,...,afc] = argmaXp,ai,...,afc/0 
subject to 



Rfe - aH - pafcaf ^ 



/c = 1,, 



(afc(s) - afc(s))^ Cfc (afc(s) - afc(s)) < 1 /c = 1, , 
Let r = 1/p. The problem (|4TI ) is equivalent to the following problem 

[f,ai,...,afc] = argmin^,ai,...,afcT 
subject to 

Rfc - cr^I afc 



K 
K. 



(41) 



.H 



>- 



k 



.K 



(42) 



^ k = l, 



.K. 



Cjt (aA;(s) - aA;(s)) 

(afc(s)-afc(s))^ 1 

This problem is once again a semi-definite programming problem that can be solved efficiently via interior 
point techniques pTJ . We can now replace the MVDR estimator by this robust version. Interestingly we 
obtain estimates of the corrected array response a(sfc). Using the model we obtain for each k 



afc(s) = ra(s). 
Hence, the self-calibration coefficients can be estimated using least squares fitting 

L 

f fc = arg min ^^ ||afc(s£) - Va.k{st)f 



7iv,7p 



(43) 



(44) 



t=i 



where F = diag{7i, ..., 7p}. Of course, when the self-calibration parameters vary slowly we can combine 
the estimation over multiple epochs. This might prove instrumental in calibration of LOFAR type arrays, 
where the calibration coefficients vary across the sky. Since the computational complexity of the self 
calibration semi-definite programming is higher than that of the MVDR dirty image, it is too complicated 
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to solve this problem for each pixel in the image. Hence it should be used similarly to the external self 
calibration cycle, where this problem is solved using a nominal source locations model. The advantage 
over ordinary self calibration, is that beyond the re-evaluation of the calibration parameters, we obtain 
better estimates of the sources powers, without significant increase in the complexity. Another interesting 
alternative proposed by the anonymous reviewer is to use the doubly constrained RCB which combines 
norm constraint like in the AAR dirty image with the robust Capon beamforming |48|. The extension to 
the moving array is done similarly to the previous problems, and will be omitted. 

IV. Statistical analysis of the 2-D MVDR DOA estimator with a moving array 

In this section we analyze the two dimensional DOA estimator based on MVDR with a moving array. 
The main motivation for this analysis is a first step in analyzing the LS-MVI, but the results have 
independent value for DOA estimation with a moving array. Interestingly, this is not a straightforward 
extension of the analysis of the 2-D MVDR estimator with a fixed array by Hawkes and Nehorai Ii49il . 
The proof is given in an appendix. 

Vaidyanathan and Buckley ||50l studied the statistical properties of 1-D MVDR estimator and a fixed 
array. They have used a single correlation matrix, and the unknown parameter was the i-th source's 
location, represented by the scalar 9i, i = 1, ...,d. Hawkes and Nehorai P9l considered the case of 2-D 
MVDR estimator and a fixed array, i.e. they have also used a single correlation matrix, but the unknown 
location was represented by the vector 6i of size (2 x 1). Our study extends these works to the case of 
a moving array. The location of the i-th source is denoted by Sj = (/j,mj), i = I, ...,d, where /j and 
rrii are the coordinates on a plane that is as an approximation of a small region in the celestial sphere 
(As described in Section JIl). R^,. is the correlation matrix at epoch k where k = l,...,K and R^ is 
the sample covariance matrix at epoch k. For simplicity we assume that for each k, R^ is based on A^ 
samples (independent of k). This is realistic in most applications. 

The 2-D MVDR spectral estimator is given by: 

t^i^ki^)^k afc(s) 
When K = I, the function degenerates to 

^^^ ~ a«(s)R-ia(s) " ^"^^^ 



Since maximizing ^ -^ is equivalent to minimizing l//(s) = a'^(s)R. ^a(s), both JSOll . ||49| 

3. (S)Jrx. 3.(S) 

minimized l//(s) for scalar and vector s, respectively. However, this approach can not be used with 

Submitted: Feb. 1st 2008, Revised: June 11th 2008. DRAFT 



18 



K covariance matrices, because of the structure of /(s) in (|45] ). In order to be able to generalize the 
analysis from a single correlation matrix to K correlation matrices, we have to directly maximize (|45]) . 
This technique significantly complicates the analysis. 

The estimator based on A^ samples per covariance matrix is denoted by (siv)j = ((^Af)i, ("iAf)j)- The 
estimation error is given by 

Asat = Sat — s. (47) 



Let s = {l,fh\ = limjv^oo ^N be the limit of the sequence of estimates. Based on general estimation 
theory the limit exists with probability one. Similarly to liSOJ , we decompose the estimation error as 



Asat = Asat + As, 



(48) 



where As = s — s is the asymptotic bias and Asat = sat — s is the finite sample error. Note that As is a 
deterministic constant, while Asat is a random variable. The variance of the finite sample error (as well 
as the variance of the estimator) decays as O (-^) and Asa? converges to with probabiUty 1. ExpUcitly 



Var[{lN)i\ = Var[{'mN)i] 



o<^'' 



This implies that the variance of the 2D MVDR estimator converges to zero asymptotically. The asymp- 
totic bias can be expUcitly computed. It is given by the following theorem: 

Theorem 4.1: Let Asj = (A/j, Arrij) be the asymptotic bias of the MVDR estimator of the i-th source's 
location, i = 1, ...,(i. 

Then, 



A/,;~ 



1^ 
1^1 



and Am,- 



1^ 
\A\ ' 



(49) 



where 



1^1 



l^il 



\A2 



dP dldm 

dldm dm? 

8/(s.,{R-^}f^J aV(s..{R-^}^lJ 
dl dldm 

8/(s.,{R-nf^J aV(s.,{R-^}f^J 
dm dm? 



aV(s.,{R-^}f^j 
aV(s.,{R-^}^^,) 

dldm 



9/(s.,{R.'nLi) 
dl 

8/(s.,{R-^}f^,) 
dm 



a/(s„{R-nf,j 

dl 



y. 47rIm{Mk,2) 
fc=i 



(M 



(50) 



kAj 
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and 



9/(s.,{Rr}f=i 



dm 






K 

E 



47r/m(Mfc,: 






fc=i 



^K 



9m^ 



yK 



k,3) 



dldm 



-7rMfc,ii?e(Mfc,6 - Mfcj)] , 

Ef=i j^WmiMk,2)Im{Mk,3) 
-Mk,iRe{Mk,8 - Mfc,9)} , 



Mfc,i = a|^R/afc 



Mfc,6 = af VfeR^ ^Vfcafc 
Mfc,2 = af UfcR^^afc Mfc,7 = af R^ ^VfcVfcafc 

Mfc,3 = a^^VfcR^^afc M^^g = af U^R^ ^V^afc 



.-1, 



Mfc,4 = a^^UfcR;, ^Ufcafc Mfc,9 = af UfcVfcR^ ^a^ 



M 



fc,5 



a^R, ^UfeUfcafc 



afe = afc(/,m) 



Ut 



exp[-27ri(uii(tfc)/ + vu{tk)m)] 

exp[-2m{u.pi{tk)l + i'pi(tfe)m)] 
uii(tfe) 



V, 







Vll{tk) 







Upl{tk) 



Vpl{tk 



(51) 



(52) 



(53) 



(54) 



(55) 



(56) 



(57) 



(58) 



The analytical expressions seem to be in a good agreement with empirical values, as demonstrated in 
section IV-BI The proof of the theorem is given in the appendix. 
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V. Simulations 

In this section we present simulation results of the LS-MVI as well as of the asymptotic bias analysis. 
In the first subsection we present deconvolution results on simulated images. We compare the performance 
of the LS-MVI algorithm and the CLEAN method. In the second subsection we study the asymptotic 
bias analysis. In all simulations we have used an East-West array of ten antenna elements, logarithmically 
spaced up to lOOOA and generated artificial sky images. We converted the image to {l,m) coordinates 
from right ascension and declination, as described in ||8l. Note that this results in {l,m) coordinates that 
are not on a rectangular grid. We produced 720 correlation matrices along 12 hours, using the model of 
|[6ll . where the averaging period for each matrix was one minute. To simplify the simulations we assumed 
perfect coherence of the sources along each integration time and that compensation for the geometric 
delay has been done. The synthesized beam used in this paper is depicted in Figure [3] at the right bottom 
of each subfigure. 

A. Comparison of LS-MVI and CLEAN 

Our first set of experiments compared the performance of the CLEAN algorithm and the LS-MVI. A 
Gaussian noise was added to the measurements. To implement the CLEAN we created a uniform grid 
and interpolated the visibilities to this grid using standard convolutional gridding. We created the classical 
dirty image using fast Fourier inversion on the rectangular grid. 

We applied the CLEAN algorithm on the created dirty image. Subsequently, we created from the 
same data the MVDR dirty image (Def. 13.11 ) and the AAR dirty image (Def. 13.2b . based on the computed 
correlation matrices. We performed deconvolution using the LS-MVI algorithm with the AAR dirty image, 
as described in Table I. Both CLEAN and LS-MVI used a loop gain 7 = 0.1 in all simulations. For each 
trial the original image, the dirty image /d(s), the MVDR dirty image /^(s) and the AAR dirty image 
I" d{s) are presented. Then the results of applying the CLEAN method and the LS-MVI algorithm. The 
latter was based on the AAR dirty image. 

The first scenario is depicted in Figure [T] The original image consists of two closely spaced sources. 
Their intensities are equal and the noise was very weak, so that the image noise was 250,000 times 
weaker than each source. It is shown that LS-MVI succeeds in separating the two sources whilst the 
CLEAN method does not. Moreover, the intensities are more accurate in the LS-MVI's image. Note that 
in this case the MVDR dirty image and the AAR dirty image are quite similar. 

The second scenario is depicted in Figure |2] The original image consists of many sources with different 
intensities Fig j2(a)[ The intensities ratio between the strongest source and the weakest source is 16dB. 
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I 




(a) True image 



(b) Dirty image 



(c) MVDR dirty image 






(d) AAR dirty image 



(e) CLEAN image 



(f) LS-MVI image 



Fig. 1. Two very strong close sources. Signal power to noise standard deviation in the image was 250,000:1. (a) The true 
image, (b) The conventional dirty image, (c) The MVDR dirty image, (d) The AAR dirty image, (e) Image after CLEAN, (f) 
Image after LS-MVI. 



The weakest source contributed to each baseline power that is equal to la of the noise standard deviation 
on the baseline. For the WSRT this is equivalent to a source of 8 mJy when using 160 MHz band and 
integrating over 256 correlation lags q As can be seen, LS-MVI performs better than CLEAN. The latter 
cannot differentiate between close sources, and the structures of its sources are less accurate. Moreover, 
LS-MVI gives better intensities estimations than CLEAN. The reason for that is the isotropic noise 
spectrum, that prevents direction dependent noise effects. Thus the sources intensities estimations are 
more accurate when based on the AAR dirty image, rather than on the conventional dirty image or on 
the MVDR dirty image. 

The third experiment included an extended source containing a central weak point source and two 
extended radio lobes. The correlation matrices were generated such that the contribution of the weak 
source to each baseline was 3(7 of the noise on the baseline. The figure presents the original image, the 
CLEAN image and the LS-MVI using the AAR dirty image. The CLEAN is presented after 100 and 120 
iterations while the LS-MVI is presented after 100 and 300 iterations respectively. Each contour is 2% 



^Based on the WSRT noise exposure calculator |http://www.astron.nI~oosterlo/expCalc.htmI| and using the fact that the 
simulated beam is of size approximately 4 pixels 
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(a) True image 



(b) Dirty image 



(c) MVDR dirty image 






(d) AAR dirty image 



(e) CLEAN image 



(f) LS-MVI image 



Fig. 2. Scattered sources. Intensities ratio between the strongest source and tiie weakest source is 16db. (a) The true image, 
(b) The conventional dirty image, (c) The MVDR dirty image, (d) The AAR dirty image, (e) Image after CLEAN, (f) Image 
after LS-MVI. 



of the dynamic range of the image. We can clearly see that the CLEAN develops a fake source a Uttle 
bit below the right lobe after 120 iteration. This causes complete divergence if the CLEAN iterations are 
continued. On the other hand, the LS-MVI with its better power estimate can continue with many more 
iterations, without affecting the image. Finally to demonstrate the improved power estimate, we have 
taken a cross section of the images near the center. Figure |4] presents the results, for the original image, 
the CLEAN after 100 iterations and the LS-MVI after 100 iterations. We can clearly see the improved 
power estimate. 

B. Bias of the 2D MVDR DOA estimator's simulation 

In this section we present a simulation of the bias of the 2D MVDR DOA estimator with a moving 
array. We have used the same array as in the previous simulations. The artificial sky image was composed 
of two point sources with different intensities (the weak source had power that was half of the power 
of the strong source). The sources contribution to each baseline was approximately lAa and 28(7 of the 
noise on the baseline. The MVDR dirty image (Def. 13.11 ) was created based on these correlation matrices. 
First we found the location of the maximum intensity in the MVDR dirty image. Then we used a 2-D 
quadratic interpolation in order to obtain a fine 2-D MVDR location estimate. We averaged the estimates 
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I 




(a) True image 



(b) CLEAN image (100 iterations) 



;^(c»i 



i 



i 



(c) CLEAN image (120 iterations) 



(d) LS-MVI image (100 iterations) 




(e) LS-MVI image (300 iterations) 

Fig. 3. Extended source. Power for weak source in the center is 3a per baseline. Contours of the image are 2% of the 
maximum. The bottom right is the synthesized beam with contours that are 20% of the maximum of the synthesized beam, (a) 
The true image, (b) Image after CLEAN - 100 iterations, (b) Image after CLEAN - 120 iterations, (c) Image after LS-MVI - 
100 iterations, (c) Image after LS-MVI - 300 iterations. 



over 100 independent trials.We compared our results to the analytical expressions given in Theorem 14. II 
This experiment was repeated 40 times with various angular separation of the sources. Figure [5] depicts 
the simulated bias (Al) against the analytical bias, as a function of the angular distance between the 
two sources in 1 coordinates. It can be seen that the analytical results are in good agreement with the 
simulation results. Moreover, as the distance between the two sources increases, the bias of the estimator 
decreases, as expected. As the distance is larger, the influence of the sidelobes of one source on the other 
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500 1 000 1 500 



Fig. 4. Cross-section through the images. 




Fig. 5. Bias of the 2-D MVDR estimator 

is smaller. The change in the bias is not monotonic in the angular separation, since the weak source 
moves through the sidelobes of the strong source. 



VI. Conclusions 

In this paper we extend the matrix formulation of [5] to non co-planar arrays and polarimetric 
measurements. Then we propose a new parametric imaging technique that improves the resolution and 
sensitivity over the classical CLEAN algorithm. The method is based on several improvements: A new 
type of dirty image, LS estimation of the powers and semi-definite constraints. We show how the technique 
can be combined into self-calibration using semi-definite programming. Our semi-definite self-calibration 
algorithm also provides a new approach to robust beamforming with a moving array, extending the 
techniques of 1401 . flU and B2l . We provide statistical analysis of the location estimator. Simulated 
examples comparing full deconvolution using LS-MVI and comparing them to the CLEAN method are 
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presented. These simulations demonstrate that the parametric approach has higher resolution, is more 
robust to source structures, and performs better in noisy situations. The great potential of the methods 
proposed in this paper is a first step, towards the development of more advanced imaging techniques, 
capable of providing higher dynamic range and interference immunity as required by the radio telescopes 
of the future. 
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VII. Appendix 
In this section we derive approximate expressions for the asymptotic bias by using a first-order Taylor 

df(s IR^^I'^ ) df(s IR^^T^ ) 

series expansion of ' gf and 'dm "" around Sj. To reduce the notational load we present 

the derivation of the single covariance matrix case. In the multiple covariance matrices case / is replaced 
by a sum and each expression can be replaced by the appropriate sum over k. The final results in theorem 
14. II includes the adjustment to multiple covariance matrices. Note that R is a fixed parameter representing 
the true covariance matrix, since we discuss only the asymptotic term. Let R^^ be fixed. The expansions 
around s is given by: 

ai ~ ai ~^ ap ' ^'* 

and 

dfi _ 

am am "*" am 

awm 
Since by definition s is the minimum of /(sj,R^^ 



ai ~ ai ^ ap • .^g^ 



df{s„R-^) ^ a/(s„R-^) aV(s„R-^) . ^^. 

am am am? * r£.iw 

+ ^^g!,'g ) -A/, + o(|(A/,|) + o(|Am.|). 



0/(s.,R-i) _ p a/(s.,R-^) _ p 

ai ~ "-'' am ~ ^■ 
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Hence, 

aV(s„R-i) A/ I aV(s„R-^) A^ ^, 9/(s..R-^) 

olom ' am^ ' am 

Proving that — ))mdi — ~ — diSm — ^^'^ using Cramer's rule, we obtain Eqs. ( |49l ) (for k=l). 

For simplicity, from now on we will denote a(sj) by a. 
Let us now derive the expressions for the derivatives of /(si,R^^). 



Therefore, 

g/(s,,R-i) _ -|[a^R-ia] 

dl [a^R-ia]2 

where 



(63) 



|[a«R-al.^R-.a+a''R-.|. (64, 



Using the fact R ^ is Hermitian and simple algebraic manipulation yields 

^R-.a + a«R-.f^2.,e(a«R-.|). 
Therefore, 



^[a^R-ia] = 2i?e(a^R-i^) . (65) 



Using Eq. ( 1651 ). Eq. (1631 ) becomes 

9/(si,R-i) ATTlm{M2) 



(66) 



9? (Mi)2 ' 

where we have used the fact that both U and R^^ are Hermitian, Mi and M2 are defined in Eq. ( [55l ) 
for k=l. Thus we get Eq. (l50l ) (for k=l). Similarly we derive Eq. (|5TI ). 
To calculate — ^-^^p (Eq. (l52l)). define (7 as 

ff(s„R-i) = 2i?e(a^R"i^). (67) 



Then, 



(68) 



5^ (a^R-ia)2 

and 

01 -2i?e(a R ^ + ^^R- ^) 
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Therefore, 

-(a^R-ia)i?e(a^R-if? + ^R-^f) • 
Using again the fact that both U and R^^ are Hermitian, we get 

^^^%^ = ^[2Im2(M2) - 7rMii?e(M4 - M5)] , (70) 

where Mi, M2, M4 and M5 are defined in Eq. dSSjl for k=l. Thus we get Eq. §^ for k=l. Eq. §3 
may be derived in the same way. 

To calculate — f)igrn — (^1- dSUl)- we use Eq. (|67] ) and obtain 

-^ '- = 2Re{--K-~'— + a''K-~'——). (71) 

om om oL oLom 

By dSHll and dlB, we get 

^%^ = (i..^w[4i^e(a^R-t)i?e(a^R-^) 
-(a^R-ia)iie(^R-it + a^R-i|;|,)]. 



(72) 



After some simple algebraic manipulations we get Eq. (1541 ) for k=l. Note that -q^^ is continuous, as 

long as a^R^^a does not vanish. Indeed a 7^ and R^^ is positive definite (since R is positive definite), 

aV(-) _ aV(-) 

dldm dmdl 



SO a R^ a 7^ 0. Therefore, q^qJ^ is continuous. Similarly, -q^^ is continuous. Hence, — -^ — — -^ 



This completes our proof. D 
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